home *** CD-ROM | disk | FTP | other *** search
/ MACD 5 / MACD 5.bin / workbench / libs / intoids.lha / Intoids 1.0 / GNU Originals / Integer.cc next >
C/C++ Source or Header  |  1995-06-16  |  47KB  |  2,283 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  16. */
  17.  
  18. /*
  19.   Some of the following algorithms are very loosely based on those from 
  20.   MIT C-Scheme bignum.c, which is
  21.       Copyright (c) 1987 Massachusetts Institute of Technology
  22.  
  23.   with other guidance from Knuth, vol. 2
  24.  
  25.   Thanks to the creators of the algorithms.
  26. */
  27.  
  28. #ifdef __GNUG__
  29. #pragma implementation
  30. #endif
  31. #include <Integer.h>
  32. #include <std.h>
  33. #include <ctype.h>
  34. #include <limits.h>
  35. #include <Obstack.h>
  36. #include <AllocRing.h>
  37. #include <new.h>
  38. #include <builtin.h>
  39. #include "Integer.hP"
  40.  
  41. #undef OK
  42.  
  43. IntRep _ZeroRep = {1, 0, 1, {0}};
  44. IntRep _OneRep = {1, 0, 1, {1}};
  45. IntRep _MinusOneRep = {1, 0, 0, {1}};
  46.  
  47.  
  48. // utilities to extract and transfer bits 
  49.  
  50. // get low bits
  51.  
  52. inline static unsigned short extract(unsigned long x)
  53. {
  54.   return x & I_MAXNUM;
  55. }
  56.  
  57. // transfer high bits to low
  58.  
  59. inline static unsigned long down(unsigned long x)
  60. {
  61.   return (x >> I_SHIFT) & I_MAXNUM;
  62. }
  63.  
  64. // transfer low bits to high
  65.  
  66. inline static unsigned long up(unsigned long x)
  67. {
  68.   return x << I_SHIFT;
  69. }
  70.  
  71. // compare two equal-length reps
  72.  
  73. inline static int docmp(const unsigned short* x, const unsigned short* y, int l)
  74. {
  75.   int diff = 0;
  76.   const unsigned short* xs = &(x[l]);
  77.   const unsigned short* ys = &(y[l]);
  78.   while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  79.   return diff;
  80. }
  81.  
  82. // figure out max length of result of +, -, etc.
  83.  
  84. inline static int calc_len(int len1, int len2, int pad)
  85. {
  86.   return (len1 >= len2)? len1 + pad : len2 + pad;
  87. }
  88.  
  89. // ensure len & sgn are correct
  90.  
  91. inline static void Icheck(IntRep* rep)
  92. {
  93.   int l = rep->len;
  94.   const unsigned short* p = &(rep->s[l]);
  95.   while (l > 0 && *--p == 0) --l;
  96.   if ((rep->len = l) == 0) rep->sgn = I_POSITIVE;
  97. }
  98.  
  99.  
  100. // zero out the end of a rep
  101.  
  102. inline static void Iclear_from(IntRep* rep, int p)
  103. {
  104.   unsigned short* cp = &(rep->s[p]);
  105.   const unsigned short* cf = &(rep->s[rep->len]);
  106.   while(cp < cf) *cp++ = 0;
  107. }
  108.  
  109. // copy parts of a rep
  110.  
  111. static inline void scpy(const unsigned short* src, unsigned short* dest,int nb)
  112. {
  113.   while (--nb >= 0) *dest++ = *src++;
  114. }
  115.  
  116. // make sure an argument is valid
  117.  
  118. static inline void nonnil(const IntRep* rep)
  119. {
  120.   if (rep == 0) 
  121.     (*lib_error_handler)("Integer", "operation on uninitialized Integer");
  122. }
  123.  
  124. // allocate a new Irep. Pad to something close to a power of two.
  125.  
  126. inline static IntRep* Inew(int newlen)
  127. {
  128.   unsigned int siz = sizeof(IntRep) + newlen * sizeof(short) + 
  129.     MALLOC_MIN_OVERHEAD;
  130.   unsigned int allocsiz = MINIntRep_SIZE;
  131.   while (allocsiz < siz) allocsiz <<= 1;  // find a power of 2
  132.   allocsiz -= MALLOC_MIN_OVERHEAD;
  133.   if (allocsiz >= MAXIntRep_SIZE * sizeof(short))
  134.     (*lib_error_handler)("Integer", "Requested length out of range");
  135.     
  136.   IntRep* rep = new (operator new (allocsiz)) IntRep;
  137.   rep->sz = (allocsiz - sizeof(IntRep) + sizeof(short)) / sizeof(short);
  138.   return rep;
  139. }
  140.  
  141. // allocate: use the bits in src if non-null, clear the rest
  142.  
  143. IntRep* Ialloc(IntRep* old, const unsigned short* src, int srclen, int newsgn,
  144.               int newlen)
  145. {
  146.   IntRep* rep;
  147.   if (old == 0 || newlen > old->sz)
  148.     rep = Inew(newlen);
  149.   else
  150.     rep = old;
  151.  
  152.   rep->len = newlen;
  153.   rep->sgn = newsgn;
  154.  
  155.   scpy(src, rep->s, srclen);
  156.   Iclear_from(rep, srclen);
  157.  
  158.   if (old != rep && old != 0 && !STATIC_IntRep(old)) delete old;
  159.   return rep;
  160. }
  161.  
  162. // allocate and clear
  163.  
  164. IntRep* Icalloc(IntRep* old, int newlen)
  165. {
  166.   IntRep* rep;
  167.   if (old == 0 || newlen > old->sz)
  168.   {
  169.     if (old != 0 && !STATIC_IntRep(old)) delete old;
  170.     rep = Inew(newlen);
  171.   }
  172.   else
  173.     rep = old;
  174.  
  175.   rep->len = newlen;
  176.   rep->sgn = I_POSITIVE;
  177.   Iclear_from(rep, 0);
  178.  
  179.   return rep;
  180. }
  181.  
  182. // reallocate
  183.  
  184. IntRep* Iresize(IntRep* old, int newlen)
  185. {
  186.   IntRep* rep;
  187.   unsigned short oldlen;
  188.   if (old == 0)
  189.   {
  190.     oldlen = 0;
  191.     rep = Inew(newlen);
  192.     rep->sgn = I_POSITIVE;
  193.   }
  194.   else 
  195.   {
  196.     oldlen = old->len;
  197.     if (newlen > old->sz)
  198.     {
  199.       rep = Inew(newlen);
  200.       scpy(old->s, rep->s, oldlen);
  201.       rep->sgn = old->sgn;
  202.       if (!STATIC_IntRep(old)) delete old;
  203.     }
  204.     else
  205.       rep = old;
  206.   }
  207.  
  208.   rep->len = newlen;
  209.   Iclear_from(rep, oldlen);
  210.  
  211.   return rep;
  212. }
  213.  
  214.  
  215. // same, for straight copy
  216.  
  217. IntRep* Icopy(IntRep* old, const IntRep* src)
  218. {
  219.   if (old == src) return old; 
  220.   IntRep* rep;
  221.   if (src == 0)
  222.   {
  223.     if (old == 0)
  224.       rep = Inew(0);
  225.     else
  226.     {
  227.       rep = old;
  228.       Iclear_from(rep, 0);
  229.     }
  230.     rep->len = 0;
  231.     rep->sgn = I_POSITIVE;
  232.   }
  233.   else 
  234.   {
  235.     int newlen = src->len;
  236.     if (old == 0 || newlen > old->sz)
  237.     {
  238.       if (old != 0 && !STATIC_IntRep(old)) delete old;
  239.       rep = Inew(newlen);
  240.     }
  241.     else
  242.       rep = old;
  243.  
  244.     rep->len = newlen;
  245.     rep->sgn = src->sgn;
  246.  
  247.     scpy(src->s, rep->s, newlen);
  248.   }
  249.  
  250.   return rep;
  251. }
  252.  
  253. // allocate & copy space for a long
  254.  
  255. IntRep* Icopy_long(IntRep* old, long x)
  256. {
  257.   int newsgn = (x >= 0);
  258.   IntRep* rep = Icopy_ulong(old, newsgn ? x : -x);
  259.   rep->sgn = newsgn;
  260.   return rep;
  261. }
  262.  
  263. IntRep* Icopy_ulong(IntRep* old, unsigned long x)
  264. {
  265.   unsigned short src[SHORT_PER_LONG];
  266.   
  267.   unsigned short srclen = 0;
  268.   while (x != 0)
  269.   {
  270.     src[srclen++] = extract(x);
  271.     x = down(x);
  272.   }
  273.  
  274.   IntRep* rep;
  275.   if (old == 0 || srclen > old->sz)
  276.   {
  277.     if (old != 0 && !STATIC_IntRep(old)) delete old;
  278.     rep = Inew(srclen);
  279.   }
  280.   else
  281.     rep = old;
  282.  
  283.   rep->len = srclen;
  284.   rep->sgn = I_POSITIVE;
  285.  
  286.   scpy(src, rep->s, srclen);
  287.  
  288.   return rep;
  289. }
  290.  
  291. // special case for zero -- it's worth it!
  292.  
  293. IntRep* Icopy_zero(IntRep* old)
  294. {
  295.   if (old == 0 || STATIC_IntRep(old))
  296.     return &_ZeroRep;
  297.  
  298.   old->len = 0;
  299.   old->sgn = I_POSITIVE;
  300.  
  301.   return old;
  302. }
  303.  
  304. // special case for 1 or -1
  305.  
  306. IntRep* Icopy_one(IntRep* old, int newsgn)
  307. {
  308.   if (old == 0 || 1 > old->sz)
  309.   {
  310.     if (old != 0 && !STATIC_IntRep(old)) delete old;
  311.     return newsgn==I_NEGATIVE ? &_MinusOneRep : &_OneRep;
  312.   }
  313.  
  314.   old->sgn = newsgn;
  315.   old->len = 1;
  316.   old->s[0] = 1;
  317.  
  318.   return old;
  319. }
  320.  
  321. // convert to a legal two's complement long if possible
  322. // if too big, return most negative/positive value
  323.  
  324. long Itolong(const IntRep* rep)
  325. { 
  326.   if ((unsigned)(rep->len) > (unsigned)(SHORT_PER_LONG))
  327.     return (rep->sgn == I_POSITIVE) ? LONG_MAX : LONG_MIN;
  328.   else if (rep->len == 0)
  329.     return 0;
  330.   else if ((unsigned)(rep->len) < (unsigned)(SHORT_PER_LONG))
  331.   {
  332.     unsigned long a = rep->s[rep->len-1];
  333.     if (SHORT_PER_LONG > 2) // normally optimized out
  334.     {
  335.       for (int i = rep->len - 2; i >= 0; --i)
  336.         a = up(a) | rep->s[i];
  337.     }
  338.     return (rep->sgn == I_POSITIVE)? a : -((long)a);
  339.   }
  340.   else 
  341.   {
  342.     unsigned long a = rep->s[SHORT_PER_LONG - 1];
  343.     if (a >= I_MINNUM)
  344.       return (rep->sgn == I_POSITIVE) ? LONG_MAX : LONG_MIN;
  345.     else
  346.     {
  347.       a = up(a) | rep->s[SHORT_PER_LONG - 2];
  348.       if (SHORT_PER_LONG > 2)
  349.       {
  350.         for (int i = SHORT_PER_LONG - 3; i >= 0; --i)
  351.           a = up(a) | rep->s[i];
  352.       }
  353.       return (rep->sgn == I_POSITIVE)? a : -((long)a);
  354.     }
  355.   }
  356. }
  357.  
  358. // test whether op long() will work.
  359. // careful about asymmetry between LONG_MIN & LONG_MAX
  360.  
  361. int Iislong(const IntRep* rep)
  362. {
  363.   unsigned int l = rep->len;
  364.   if (l < SHORT_PER_LONG)
  365.     return 1;
  366.   else if (l > SHORT_PER_LONG)
  367.     return 0;
  368.   else if ((unsigned)(rep->s[SHORT_PER_LONG - 1]) < (unsigned)(I_MINNUM))
  369.     return 1;
  370.   else if (rep->sgn == I_NEGATIVE && rep->s[SHORT_PER_LONG - 1] == I_MINNUM)
  371.   {
  372.     for (unsigned int i = 0; i < SHORT_PER_LONG - 1; ++i)
  373.       if (rep->s[i] != 0)
  374.         return 0;
  375.     return 1;
  376.   }
  377.   else
  378.     return 0;
  379. }
  380.  
  381. // comparison functions
  382.   
  383. int compare(const IntRep* x, const IntRep* y)
  384. {
  385.   int diff  = x->sgn - y->sgn;
  386.   if (diff == 0)
  387.   {
  388.     diff = x->len - y->len;
  389.     if (diff == 0)
  390.       diff = docmp(x->s, y->s, x->len);
  391.     if (x->sgn == I_NEGATIVE)
  392.       diff = -diff;
  393.   }
  394.   return diff;
  395. }
  396.  
  397. int ucompare(const IntRep* x, const IntRep* y)
  398. {
  399.   int diff = x->len - y->len;
  400.   if (diff == 0)
  401.   {
  402.     int l = x->len;
  403.     const unsigned short* xs = &(x->s[l]);
  404.     const unsigned short* ys = &(y->s[l]);
  405.     while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  406.   }
  407.   return diff;
  408. }
  409.  
  410. int compare(const IntRep* x, long  y)
  411. {
  412.   int xl = x->len;
  413.   int xsgn = x->sgn;
  414.   if (y == 0)
  415.   {
  416.     if (xl == 0)
  417.       return 0;
  418.     else if (xsgn == I_NEGATIVE)
  419.       return -1;
  420.     else
  421.       return 1;
  422.   }
  423.   else
  424.   {
  425.     int ysgn = y >= 0;
  426.     unsigned long uy = (ysgn)? y : -y;
  427.     int diff = xsgn - ysgn;
  428.     if (diff == 0)
  429.     {
  430.       diff = xl - SHORT_PER_LONG;
  431.       if (diff <= 0)
  432.       {
  433.         unsigned short tmp[SHORT_PER_LONG];
  434.         int yl = 0;
  435.         while (uy != 0)
  436.         {
  437.           tmp[yl++] = extract(uy);
  438.           uy = down(uy);
  439.         }
  440.         diff = xl - yl;
  441.         if (diff == 0)
  442.           diff = docmp(x->s, tmp, xl);
  443.       }
  444.       if (xsgn == I_NEGATIVE)
  445.     diff = -diff;
  446.     }
  447.     return diff;
  448.   }
  449. }
  450.  
  451. int ucompare(const IntRep* x, long  y)
  452. {
  453.   int xl = x->len;
  454.   if (y == 0)
  455.     return xl;
  456.   else
  457.   {
  458.     unsigned long uy = (y >= 0)? y : -y;
  459.     int diff = xl - SHORT_PER_LONG;
  460.     if (diff <= 0)
  461.     {
  462.       unsigned short tmp[SHORT_PER_LONG];
  463.       int yl = 0;
  464.       while (uy != 0)
  465.       {
  466.         tmp[yl++] = extract(uy);
  467.         uy = down(uy);
  468.       }
  469.       diff = xl - yl;
  470.       if (diff == 0)
  471.         diff = docmp(x->s, tmp, xl);
  472.     }
  473.     return diff;
  474.   }
  475. }
  476.  
  477.  
  478.  
  479. // arithmetic functions
  480.  
  481. IntRep* add(const IntRep* x, int negatex, 
  482.             const IntRep* y, int negatey, IntRep* r)
  483. {
  484.   nonnil(x);
  485.   nonnil(y);
  486.  
  487.   int xl = x->len;
  488.   int yl = y->len;
  489.  
  490.   int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  491.   int ysgn = (negatey && yl != 0) ? !y->sgn : y->sgn;
  492.  
  493.   int xrsame = x == r;
  494.   int yrsame = y == r;
  495.  
  496.   if (yl == 0)
  497.     r = Ialloc(r, x->s, xl, xsgn, xl);
  498.   else if (xl == 0)
  499.     r = Ialloc(r, y->s, yl, ysgn, yl);
  500.   else if (xsgn == ysgn)
  501.   {
  502.     if (xrsame || yrsame)
  503.       r = Iresize(r, calc_len(xl, yl, 1));
  504.     else
  505.       r = Icalloc(r, calc_len(xl, yl, 1));
  506.     r->sgn = xsgn;
  507.     unsigned short* rs = r->s;
  508.     const unsigned short* as;
  509.     const unsigned short* bs;
  510.     const unsigned short* topa;
  511.     const unsigned short* topb;
  512.     if (xl >= yl)
  513.     {
  514.       as =  (xrsame)? r->s : x->s;
  515.       topa = &(as[xl]);
  516.       bs =  (yrsame)? r->s : y->s;
  517.       topb = &(bs[yl]);
  518.     }
  519.     else
  520.     {
  521.       bs =  (xrsame)? r->s : x->s;
  522.       topb = &(bs[xl]);
  523.       as =  (yrsame)? r->s : y->s;
  524.       topa = &(as[yl]);
  525.     }
  526.     unsigned long sum = 0;
  527.     while (bs < topb)
  528.     {
  529.       sum += (unsigned long)(*as++) + (unsigned long)(*bs++);
  530.       *rs++ = extract(sum);
  531.       sum = down(sum);
  532.     }
  533.     while (sum != 0 && as < topa)
  534.     {
  535.       sum += (unsigned long)(*as++);
  536.       *rs++ = extract(sum);
  537.       sum = down(sum);
  538.     }
  539.     if (sum != 0)
  540.       *rs = extract(sum);
  541.     else if (rs != as)
  542.       while (as < topa)
  543.         *rs++ = *as++;
  544.   }
  545.   else
  546.   {
  547.     int comp = ucompare(x, y);
  548.     if (comp == 0)
  549.       r = Icopy_zero(r);
  550.     else
  551.     {
  552.       if (xrsame || yrsame)
  553.         r = Iresize(r, calc_len(xl, yl, 0));
  554.       else
  555.         r = Icalloc(r, calc_len(xl, yl, 0));
  556.       unsigned short* rs = r->s;
  557.       const unsigned short* as;
  558.       const unsigned short* bs;
  559.       const unsigned short* topa;
  560.       const unsigned short* topb;
  561.       if (comp > 0)
  562.       {
  563.         as =  (xrsame)? r->s : x->s;
  564.         topa = &(as[xl]);
  565.         bs =  (yrsame)? r->s : y->s;
  566.         topb = &(bs[yl]);
  567.         r->sgn = xsgn;
  568.       }
  569.       else
  570.       {
  571.         bs =  (xrsame)? r->s : x->s;
  572.         topb = &(bs[xl]);
  573.         as =  (yrsame)? r->s : y->s;
  574.         topa = &(as[yl]);
  575.         r->sgn = ysgn;
  576.       }
  577.       unsigned long hi = 1;
  578.       while (bs < topb)
  579.       {
  580.         hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
  581.         *rs++ = extract(hi);
  582.         hi = down(hi);
  583.       }
  584.       while (hi == 0 && as < topa)
  585.       {
  586.         hi = (unsigned long)(*as++) + I_MAXNUM;
  587.         *rs++ = extract(hi);
  588.         hi = down(hi);
  589.       }
  590.       if (rs != as)
  591.         while (as < topa)
  592.           *rs++ = *as++;
  593.     }
  594.   }
  595.   Icheck(r);
  596.   return r;
  597. }
  598.  
  599.  
  600. IntRep* add(const IntRep* x, int negatex, long y, IntRep* r)
  601. {
  602.   nonnil(x);
  603.   int xl = x->len;
  604.   int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  605.   int xrsame = x == r;
  606.  
  607.   int ysgn = (y >= 0);
  608.   unsigned long uy = (ysgn)? y : -y;
  609.  
  610.   if (y == 0)
  611.     r = Ialloc(r, x->s, xl, xsgn, xl);
  612.   else if (xl == 0)
  613.     r = Icopy_long(r, y);
  614.   else if (xsgn == ysgn)
  615.   {
  616.     if (xrsame)
  617.       r = Iresize(r, calc_len(xl, SHORT_PER_LONG, 1));
  618.     else
  619.       r = Icalloc(r, calc_len(xl, SHORT_PER_LONG, 1));
  620.     r->sgn = xsgn;
  621.     unsigned short* rs = r->s;
  622.     const unsigned short* as =  (xrsame)? r->s : x->s;
  623.     const unsigned short* topa = &(as[xl]);
  624.     unsigned long sum = 0;
  625.     while (as < topa && uy != 0)
  626.     {
  627.       unsigned long u = extract(uy);
  628.       uy = down(uy);
  629.       sum += (unsigned long)(*as++) + u;
  630.       *rs++ = extract(sum);
  631.       sum = down(sum);
  632.     }
  633.     while (sum != 0 && as < topa)
  634.     {
  635.       sum += (unsigned long)(*as++);
  636.       *rs++ = extract(sum);
  637.       sum = down(sum);
  638.     }
  639.     if (sum != 0)
  640.       *rs = extract(sum);
  641.     else if (rs != as)
  642.       while (as < topa)
  643.         *rs++ = *as++;
  644.   }
  645.   else
  646.   {
  647.     unsigned short tmp[SHORT_PER_LONG];
  648.     int yl = 0;
  649.     while (uy != 0)
  650.     {
  651.       tmp[yl++] = extract(uy);
  652.       uy = down(uy);
  653.     }
  654.     int comp = xl - yl;
  655.     if (comp == 0)
  656.       comp = docmp(x->s, tmp, yl);
  657.     if (comp == 0)
  658.       r = Icopy_zero(r);
  659.     else
  660.     {
  661.       if (xrsame)
  662.         r = Iresize(r, calc_len(xl, yl, 0));
  663.       else
  664.         r = Icalloc(r, calc_len(xl, yl, 0));
  665.       unsigned short* rs = r->s;
  666.       const unsigned short* as;
  667.       const unsigned short* bs;
  668.       const unsigned short* topa;
  669.       const unsigned short* topb;
  670.       if (comp > 0)
  671.       {
  672.         as =  (xrsame)? r->s : x->s;
  673.         topa = &(as[xl]);
  674.         bs =  tmp;
  675.         topb = &(bs[yl]);
  676.         r->sgn = xsgn;
  677.       }
  678.       else
  679.       {
  680.         bs =  (xrsame)? r->s : x->s;
  681.         topb = &(bs[xl]);
  682.         as =  tmp;
  683.         topa = &(as[yl]);
  684.         r->sgn = ysgn;
  685.       }
  686.       unsigned long hi = 1;
  687.       while (bs < topb)
  688.       {
  689.         hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
  690.         *rs++ = extract(hi);
  691.         hi = down(hi);
  692.       }
  693.       while (hi == 0 && as < topa)
  694.       {
  695.         hi = (unsigned long)(*as++) + I_MAXNUM;
  696.         *rs++ = extract(hi);
  697.         hi = down(hi);
  698.       }
  699.       if (rs != as)
  700.         while (as < topa)
  701.           *rs++ = *as++;
  702.     }
  703.   }
  704.   Icheck(r);
  705.   return r;
  706. }
  707.  
  708.  
  709. IntRep* multiply(const IntRep* x, const IntRep* y, IntRep* r)
  710. {
  711.   nonnil(x);
  712.   nonnil(y);
  713.   int xl = x->len;
  714.   int yl = y->len;
  715.   int rl = xl + yl;
  716.   int rsgn = x->sgn == y->sgn;
  717.   int xrsame = x == r;
  718.   int yrsame = y == r;
  719.   int xysame = x == y;
  720.   
  721.   if (xl == 0 || yl == 0)
  722.     r = Icopy_zero(r);
  723.   else if (xl == 1 && x->s[0] == 1)
  724.     r = Icopy(r, y);
  725.   else if (yl == 1 && y->s[0] == 1)
  726.     r = Icopy(r, x);
  727.   else if (!(xysame && xrsame))
  728.   {
  729.     if (xrsame || yrsame)
  730.       r = Iresize(r, rl);
  731.     else
  732.       r = Icalloc(r, rl);
  733.     unsigned short* rs = r->s;
  734.     unsigned short* topr = &(rs[rl]);
  735.  
  736.     // use best inner/outer loop params given constraints
  737.     unsigned short* currentr;
  738.     const unsigned short* bota;
  739.     const unsigned short* as;
  740.     const unsigned short* botb;
  741.     const unsigned short* topb;
  742.     if (xrsame)                 
  743.     { 
  744.       currentr = &(rs[xl-1]);
  745.       bota = rs;
  746.       as = currentr;
  747.       botb = y->s;
  748.       topb = &(botb[yl]);
  749.     }
  750.     else if (yrsame)
  751.     {
  752.       currentr = &(rs[yl-1]);
  753.       bota = rs;
  754.       as = currentr;
  755.       botb = x->s;
  756.       topb = &(botb[xl]);
  757.     }
  758.     else if (xl <= yl)
  759.     {
  760.       currentr = &(rs[xl-1]);
  761.       bota = x->s;
  762.       as = &(bota[xl-1]);
  763.       botb = y->s;
  764.       topb = &(botb[yl]);
  765.     }
  766.     else
  767.     {
  768.       currentr = &(rs[yl-1]);
  769.       bota = y->s;
  770.       as = &(bota[yl-1]);
  771.       botb = x->s;
  772.       topb = &(botb[xl]);
  773.     }
  774.  
  775.     while (as >= bota)
  776.     {
  777.       unsigned long ai = (unsigned long)(*as--);
  778.       unsigned short* rs = currentr--;
  779.       *rs = 0;
  780.       if (ai != 0)
  781.       {
  782.         unsigned long sum = 0;
  783.         const unsigned short* bs = botb;
  784.         while (bs < topb)
  785.         {
  786.           sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
  787.           *rs++ = extract(sum);
  788.           sum = down(sum);
  789.         }
  790.         while (sum != 0 && rs < topr)
  791.         {
  792.           sum += (unsigned long)(*rs);
  793.           *rs++ = extract(sum);
  794.           sum = down(sum);
  795.         }
  796.       }
  797.     }
  798.   }
  799.   else                          // x, y, and r same; compute over diagonals
  800.   {
  801.     r = Iresize(r, rl);
  802.     unsigned short* botr = r->s;
  803.     unsigned short* topr = &(botr[rl]);
  804.     unsigned short* rs =   &(botr[rl - 2]);
  805.  
  806.     const unsigned short* bota = (xrsame)? botr : x->s;
  807.     const unsigned short* loa =  &(bota[xl - 1]);
  808.     const unsigned short* hia =  loa;
  809.  
  810.     for (; rs >= botr; --rs)
  811.     {
  812.       const unsigned short* h = hia;
  813.       const unsigned short* l = loa;
  814.       unsigned long prod = (unsigned long)(*h) * (unsigned long)(*l);
  815.       *rs = 0;
  816.  
  817.       for(;;)
  818.       {
  819.         unsigned short* rt = rs;
  820.         unsigned long sum = prod + (unsigned long)(*rt);
  821.         *rt++ = extract(sum);
  822.         sum = down(sum);
  823.         while (sum != 0 && rt < topr)
  824.         {
  825.           sum += (unsigned long)(*rt);
  826.           *rt++ = extract(sum);
  827.           sum = down(sum);
  828.         }
  829.         if (h > l)
  830.         {
  831.           rt = rs;
  832.           sum = prod + (unsigned long)(*rt);
  833.           *rt++ = extract(sum);
  834.           sum = down(sum);
  835.           while (sum != 0 && rt < topr)
  836.           {
  837.             sum += (unsigned long)(*rt);
  838.             *rt++ = extract(sum);
  839.             sum = down(sum);
  840.           }
  841.           if (--h >= ++l)
  842.             prod = (unsigned long)(*h) * (unsigned long)(*l);
  843.           else
  844.             break;
  845.         }
  846.         else
  847.           break;
  848.       }
  849.       if (loa > bota)
  850.         --loa;
  851.       else
  852.         --hia;
  853.     }
  854.   }
  855.   r->sgn = rsgn;
  856.   Icheck(r);
  857.   return r;
  858. }
  859.  
  860.  
  861. IntRep* multiply(const IntRep* x, long y, IntRep* r)
  862. {
  863.   nonnil(x);
  864.   int xl = x->len;
  865.     
  866.   if (xl == 0 || y == 0)
  867.     r = Icopy_zero(r);
  868.   else if (y == 1)
  869.     r = Icopy(r, x);
  870.   else
  871.   {
  872.     int ysgn = y >= 0;
  873.     int rsgn = x->sgn == ysgn;
  874.     unsigned long uy = (ysgn)? y : -y;
  875.     unsigned short tmp[SHORT_PER_LONG];
  876.     int yl = 0;
  877.     while (uy != 0)
  878.     {
  879.       tmp[yl++] = extract(uy);
  880.       uy = down(uy);
  881.     }
  882.  
  883.     int rl = xl + yl;
  884.     int xrsame = x == r;
  885.     if (xrsame)
  886.       r = Iresize(r, rl);
  887.     else
  888.       r = Icalloc(r, rl);
  889.  
  890.     unsigned short* rs = r->s;
  891.     unsigned short* topr = &(rs[rl]);
  892.     unsigned short* currentr;
  893.     const unsigned short* bota;
  894.     const unsigned short* as;
  895.     const unsigned short* botb;
  896.     const unsigned short* topb;
  897.  
  898.     if (xrsame)
  899.     { 
  900.       currentr = &(rs[xl-1]);
  901.       bota = rs;
  902.       as = currentr;
  903.       botb = tmp;
  904.       topb = &(botb[yl]);
  905.     }
  906.     else if (xl <= yl)
  907.     {
  908.       currentr = &(rs[xl-1]);
  909.       bota = x->s;
  910.       as = &(bota[xl-1]);
  911.       botb = tmp;
  912.       topb = &(botb[yl]);
  913.     }
  914.     else
  915.     {
  916.       currentr = &(rs[yl-1]);
  917.       bota = tmp;
  918.       as = &(bota[yl-1]);
  919.       botb = x->s;
  920.       topb = &(botb[xl]);
  921.     }
  922.  
  923.     while (as >= bota)
  924.     {
  925.       unsigned long ai = (unsigned long)(*as--);
  926.       unsigned short* rs = currentr--;
  927.       *rs = 0;
  928.       if (ai != 0)
  929.       {
  930.         unsigned long sum = 0;
  931.         const unsigned short* bs = botb;
  932.         while (bs < topb)
  933.         {
  934.           sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
  935.           *rs++ = extract(sum);
  936.           sum = down(sum);
  937.         }
  938.         while (sum != 0 && rs < topr)
  939.         {
  940.           sum += (unsigned long)(*rs);
  941.           *rs++ = extract(sum);
  942.           sum = down(sum);
  943.         }
  944.       }
  945.     }
  946.     r->sgn = rsgn;
  947.   }
  948.   Icheck(r);
  949.   return r;
  950. }
  951.  
  952.  
  953. // main division routine
  954.  
  955. static void do_divide(unsigned short* rs,
  956.                       const unsigned short* ys, int yl,
  957.                       unsigned short* qs, int ql)
  958. {
  959.   const unsigned short* topy = &(ys[yl]);
  960.   unsigned short d1 = ys[yl - 1];
  961.   unsigned short d2 = ys[yl - 2];
  962.  
  963.   int l = ql - 1;
  964.   int i = l + yl;
  965.   
  966.   for (; l >= 0; --l, --i)
  967.   {
  968.     unsigned short qhat;       // guess q
  969.     if (d1 == rs[i])
  970.       qhat = I_MAXNUM;
  971.     else
  972.     {
  973.       unsigned long lr = up((unsigned long)rs[i]) | rs[i-1];
  974.       qhat = lr / d1;
  975.     }
  976.  
  977.     for(;;)     // adjust q, use docmp to avoid overflow problems
  978.     {
  979.       unsigned short ts[3];
  980.       unsigned long prod = (unsigned long)d2 * (unsigned long)qhat;
  981.       ts[0] = extract(prod);
  982.       prod = down(prod) + (unsigned long)d1 * (unsigned long)qhat;
  983.       ts[1] = extract(prod);
  984.       ts[2] = extract(down(prod));
  985.       if (docmp(ts, &(rs[i-2]), 3) > 0)
  986.         --qhat;
  987.       else
  988.         break;
  989.     };
  990.     
  991.     // multiply & subtract
  992.     
  993.     const unsigned short* yt = ys;
  994.     unsigned short* rt = &(rs[l]);
  995.     unsigned long prod = 0;
  996.     unsigned long hi = 1;
  997.     while (yt < topy)
  998.     {
  999.       prod = (unsigned long)qhat * (unsigned long)(*yt++) + down(prod);
  1000.       hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(extract(prod));
  1001.       *rt++ = extract(hi);
  1002.       hi = down(hi);
  1003.     }
  1004.     hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(down(prod));
  1005.     *rt = extract(hi);
  1006.     hi = down(hi);
  1007.     
  1008.     // off-by-one, add back
  1009.     
  1010.     if (hi == 0)
  1011.     {
  1012.       --qhat;
  1013.       yt = ys;
  1014.       rt = &(rs[l]);
  1015.       hi = 0;
  1016.       while (yt < topy)
  1017.       {
  1018.         hi = (unsigned long)(*rt) + (unsigned long)(*yt++) + down(hi);
  1019.         *rt++ = extract(hi);
  1020.       }
  1021.       *rt = 0;
  1022.     }
  1023.     if (qs != 0)
  1024.       qs[l] = qhat;
  1025.   }
  1026. }
  1027.  
  1028. // divide by single digit, return remainder
  1029. // if q != 0, then keep the result in q, else just compute rem
  1030.  
  1031. static int unscale(const unsigned short* x, int xl, unsigned short y,
  1032.                    unsigned short* q)
  1033. {
  1034.   if (xl == 0 || y == 1)
  1035.     return 0;
  1036.   else if (q != 0)
  1037.   {
  1038.     unsigned short* botq = q;
  1039.     unsigned short* qs = &(botq[xl - 1]);
  1040.     const unsigned short* xs = &(x[xl - 1]);
  1041.     unsigned long rem = 0;
  1042.     while (qs >= botq)
  1043.     {
  1044.       rem = up(rem) | *xs--;
  1045.       unsigned long u = rem / y;
  1046.       *qs-- = extract(u);
  1047.       rem -= u * y;
  1048.     }
  1049.     int r = extract(rem);
  1050.     return r;
  1051.   }
  1052.   else                          // same loop, a bit faster if just need rem
  1053.   {
  1054.     const unsigned short* botx = x;
  1055.     const unsigned short* xs = &(botx[xl - 1]);
  1056.     unsigned long rem = 0;
  1057.     while (xs >= botx)
  1058.     {
  1059.       rem = up(rem) | *xs--;
  1060.       unsigned long u = rem / y;
  1061.       rem -= u * y;
  1062.     }
  1063.     int r = extract(rem);
  1064.     return r;
  1065.   }
  1066. }
  1067.  
  1068.  
  1069. IntRep* div(const IntRep* x, const IntRep* y, IntRep* q)
  1070. {
  1071.   nonnil(x);
  1072.   nonnil(y);
  1073.   int xl = x->len;
  1074.   int yl = y->len;
  1075.   if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1076.  
  1077.   int comp = ucompare(x, y);
  1078.   int xsgn = x->sgn;
  1079.   int ysgn = y->sgn;
  1080.  
  1081.   int samesign = xsgn == ysgn;
  1082.  
  1083.   if (comp < 0)
  1084.     q = Icopy_zero(q);
  1085.   else if (comp == 0)
  1086.     q = Icopy_one(q, samesign);
  1087.   else if (yl == 1)
  1088.   {
  1089.     q = Icopy(q, x);
  1090.     unscale(q->s, q->len, y->s[0], q->s);
  1091.   }
  1092.   else
  1093.   {
  1094.     IntRep* yy = 0;
  1095.     IntRep* r  = 0;
  1096.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1097.     if (prescale != 1 || y == q)
  1098.     {
  1099.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1100.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1101.     }
  1102.     else
  1103.     {
  1104.       yy = (IntRep*)y;
  1105.       r = Icalloc(r, xl + 1);
  1106.       scpy(x->s, r->s, xl);
  1107.     }
  1108.  
  1109.     int ql = xl - yl + 1;
  1110.       
  1111.     q = Icalloc(q, ql);
  1112.     do_divide(r->s, yy->s, yl, q->s, ql);
  1113.  
  1114.     if (yy != y && !STATIC_IntRep(yy)) delete yy;
  1115.     if (!STATIC_IntRep(r)) delete r;
  1116.   }
  1117.   q->sgn = samesign;
  1118.   Icheck(q);
  1119.   return q;
  1120. }
  1121.  
  1122. IntRep* div(const IntRep* x, long y, IntRep* q)
  1123. {
  1124.   nonnil(x);
  1125.   int xl = x->len;
  1126.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1127.  
  1128.   unsigned short ys[SHORT_PER_LONG];
  1129.   unsigned long u;
  1130.   int ysgn = y >= 0;
  1131.   if (ysgn)
  1132.     u = y;
  1133.   else
  1134.     u = -y;
  1135.   int yl = 0;
  1136.   while (u != 0)
  1137.   {
  1138.     ys[yl++] = extract(u);
  1139.     u = down(u);
  1140.   }
  1141.  
  1142.   int comp = xl - yl;
  1143.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1144.  
  1145.   int xsgn = x->sgn;
  1146.   int samesign = xsgn == ysgn;
  1147.  
  1148.   if (comp < 0)
  1149.     q = Icopy_zero(q);
  1150.   else if (comp == 0)
  1151.   {
  1152.     q = Icopy_one(q, samesign);
  1153.   }
  1154.   else if (yl == 1)
  1155.   {
  1156.     q = Icopy(q, x);
  1157.     unscale(q->s, q->len, ys[0], q->s);
  1158.   }
  1159.   else
  1160.   {
  1161.     IntRep* r  = 0;
  1162.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1163.     if (prescale != 1)
  1164.     {
  1165.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1166.       ys[0] = extract(prod);
  1167.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1168.       ys[1] = extract(prod);
  1169.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1170.     }
  1171.     else
  1172.     {
  1173.       r = Icalloc(r, xl + 1);
  1174.       scpy(x->s, r->s, xl);
  1175.     }
  1176.  
  1177.     int ql = xl - yl + 1;
  1178.       
  1179.     q = Icalloc(q, ql);
  1180.     do_divide(r->s, ys, yl, q->s, ql);
  1181.  
  1182.     if (!STATIC_IntRep(r)) delete r;
  1183.   }
  1184.   q->sgn = samesign;
  1185.   Icheck(q);
  1186.   return q;
  1187. }
  1188.  
  1189.  
  1190. void divide(const Integer& Ix, long y, Integer& Iq, long& rem)
  1191. {
  1192.   const IntRep* x = Ix.rep;
  1193.   nonnil(x);
  1194.   IntRep* q = Iq.rep;
  1195.   int xl = x->len;
  1196.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1197.  
  1198.   unsigned short ys[SHORT_PER_LONG];
  1199.   unsigned long u;
  1200.   int ysgn = y >= 0;
  1201.   if (ysgn)
  1202.     u = y;
  1203.   else
  1204.     u = -y;
  1205.   int yl = 0;
  1206.   while (u != 0)
  1207.   {
  1208.     ys[yl++] = extract(u);
  1209.     u = down(u);
  1210.   }
  1211.  
  1212.   int comp = xl - yl;
  1213.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1214.  
  1215.   int xsgn = x->sgn;
  1216.   int samesign = xsgn == ysgn;
  1217.  
  1218.   if (comp < 0)
  1219.   {
  1220.     rem = Itolong(x);
  1221.     q = Icopy_zero(q);
  1222.   }
  1223.   else if (comp == 0)
  1224.   {
  1225.     q = Icopy_one(q, samesign);
  1226.     rem = 0;
  1227.   }
  1228.   else if (yl == 1)
  1229.   {
  1230.     q = Icopy(q, x);
  1231.     rem = unscale(q->s, q->len, ys[0], q->s);
  1232.   }
  1233.   else
  1234.   {
  1235.     IntRep* r  = 0;
  1236.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1237.     if (prescale != 1)
  1238.     {
  1239.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1240.       ys[0] = extract(prod);
  1241.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1242.       ys[1] = extract(prod);
  1243.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1244.     }
  1245.     else
  1246.     {
  1247.       r = Icalloc(r, xl + 1);
  1248.       scpy(x->s, r->s, xl);
  1249.     }
  1250.  
  1251.     int ql = xl - yl + 1;
  1252.       
  1253.     q = Icalloc(q, ql);
  1254.     
  1255.     do_divide(r->s, ys, yl, q->s, ql);
  1256.  
  1257.     if (prescale != 1)
  1258.     {
  1259.       Icheck(r);
  1260.       unscale(r->s, r->len, prescale, r->s);
  1261.     }
  1262.     Icheck(r);
  1263.     rem = Itolong(r);
  1264.     if (!STATIC_IntRep(r)) delete r;
  1265.   }
  1266.   rem = abs(rem);
  1267.   if (xsgn == I_NEGATIVE) rem = -rem;
  1268.   q->sgn = samesign;
  1269.   Icheck(q);
  1270.   Iq.rep = q;
  1271. }
  1272.  
  1273.  
  1274. void divide(const Integer& Ix, const Integer& Iy, Integer& Iq, Integer& Ir)
  1275. {
  1276.   const IntRep* x = Ix.rep;
  1277.   nonnil(x);
  1278.   const IntRep* y = Iy.rep;
  1279.   nonnil(y);
  1280.   IntRep* q = Iq.rep;
  1281.   IntRep* r = Ir.rep;
  1282.  
  1283.   int xl = x->len;
  1284.   int yl = y->len;
  1285.   if (yl == 0)
  1286.     (*lib_error_handler)("Integer", "attempted division by zero");
  1287.  
  1288.   int comp = ucompare(x, y);
  1289.   int xsgn = x->sgn;
  1290.   int ysgn = y->sgn;
  1291.  
  1292.   int samesign = xsgn == ysgn;
  1293.  
  1294.   if (comp < 0)
  1295.   {
  1296.     q = Icopy_zero(q);
  1297.     r = Icopy(r, x);
  1298.   }
  1299.   else if (comp == 0)
  1300.   {
  1301.     q = Icopy_one(q, samesign);
  1302.     r = Icopy_zero(r);
  1303.   }
  1304.   else if (yl == 1)
  1305.   {
  1306.     q = Icopy(q, x);
  1307.     int rem =  unscale(q->s, q->len, y->s[0], q->s);
  1308.     r = Icopy_long(r, rem);
  1309.     if (rem != 0)
  1310.       r->sgn = xsgn;
  1311.   }
  1312.   else
  1313.   {
  1314.     IntRep* yy = 0;
  1315.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1316.     if (prescale != 1 || y == q || y == r)
  1317.     {
  1318.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1319.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1320.     }
  1321.     else
  1322.     {
  1323.       yy = (IntRep*)y;
  1324.       r = Icalloc(r, xl + 1);
  1325.       scpy(x->s, r->s, xl);
  1326.     }
  1327.  
  1328.     int ql = xl - yl + 1;
  1329.       
  1330.     q = Icalloc(q, ql);
  1331.     do_divide(r->s, yy->s, yl, q->s, ql);
  1332.  
  1333.     if (yy != y && !STATIC_IntRep(yy)) delete yy;
  1334.     if (prescale != 1)
  1335.     {
  1336.       Icheck(r);
  1337.       unscale(r->s, r->len, prescale, r->s);
  1338.     }
  1339.   }
  1340.   q->sgn = samesign;
  1341.   Icheck(q);
  1342.   Iq.rep = q;
  1343.   Icheck(r);
  1344.   Ir.rep = r;
  1345. }
  1346.  
  1347. IntRep* mod(const IntRep* x, const IntRep* y, IntRep* r)
  1348. {
  1349.   nonnil(x);
  1350.   nonnil(y);
  1351.   int xl = x->len;
  1352.   int yl = y->len;
  1353.   if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1354.  
  1355.   int comp = ucompare(x, y);
  1356.   int xsgn = x->sgn;
  1357.  
  1358.   if (comp < 0)
  1359.     r = Icopy(r, x);
  1360.   else if (comp == 0)
  1361.     r = Icopy_zero(r);
  1362.   else if (yl == 1)
  1363.   {
  1364.     int rem = unscale(x->s, xl, y->s[0], 0);
  1365.     r = Icopy_long(r, rem);
  1366.     if (rem != 0)
  1367.       r->sgn = xsgn;
  1368.   }
  1369.   else
  1370.   {
  1371.     IntRep* yy = 0;
  1372.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1373.     if (prescale != 1 || y == r)
  1374.     {
  1375.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1376.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1377.     }
  1378.     else
  1379.     {
  1380.       yy = (IntRep*)y;
  1381.       r = Icalloc(r, xl + 1);
  1382.       scpy(x->s, r->s, xl);
  1383.     }
  1384.       
  1385.     do_divide(r->s, yy->s, yl, 0, xl - yl + 1);
  1386.  
  1387.     if (yy != y && !STATIC_IntRep(yy)) delete yy;
  1388.  
  1389.     if (prescale != 1)
  1390.     {
  1391.       Icheck(r);
  1392.       unscale(r->s, r->len, prescale, r->s);
  1393.     }
  1394.   }
  1395.   Icheck(r);
  1396.   return r;
  1397. }
  1398.  
  1399. IntRep* mod(const IntRep* x, long y, IntRep* r)
  1400. {
  1401.   nonnil(x);
  1402.   int xl = x->len;
  1403.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1404.  
  1405.   unsigned short ys[SHORT_PER_LONG];
  1406.   unsigned long u;
  1407.   int ysgn = y >= 0;
  1408.   if (ysgn)
  1409.     u = y;
  1410.   else
  1411.     u = -y;
  1412.   int yl = 0;
  1413.   while (u != 0)
  1414.   {
  1415.     ys[yl++] = extract(u);
  1416.     u = down(u);
  1417.   }
  1418.  
  1419.   int comp = xl - yl;
  1420.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1421.  
  1422.   int xsgn = x->sgn;
  1423.  
  1424.   if (comp < 0)
  1425.     r = Icopy(r, x);
  1426.   else if (comp == 0)
  1427.     r = Icopy_zero(r);
  1428.   else if (yl == 1)
  1429.   {
  1430.     int rem = unscale(x->s, xl, ys[0], 0);
  1431.     r = Icopy_long(r, rem);
  1432.     if (rem != 0)
  1433.       r->sgn = xsgn;
  1434.   }
  1435.   else
  1436.   {
  1437.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1438.     if (prescale != 1)
  1439.     {
  1440.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1441.       ys[0] = extract(prod);
  1442.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1443.       ys[1] = extract(prod);
  1444.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1445.     }
  1446.     else
  1447.     {
  1448.       r = Icalloc(r, xl + 1);
  1449.       scpy(x->s, r->s, xl);
  1450.     }
  1451.       
  1452.     do_divide(r->s, ys, yl, 0, xl - yl + 1);
  1453.  
  1454.     if (prescale != 1)
  1455.     {
  1456.       Icheck(r);
  1457.       unscale(r->s, r->len, prescale, r->s);
  1458.     }
  1459.   }
  1460.   Icheck(r);
  1461.   return r;
  1462. }
  1463.  
  1464. IntRep* lshift(const IntRep* x, long y, IntRep* r)
  1465. {
  1466.   nonnil(x);
  1467.   int xl = x->len;
  1468.   if (xl == 0 || y == 0)
  1469.   {
  1470.     r = Icopy(r, x);
  1471.     return r;
  1472.   }
  1473.  
  1474.   int xrsame = x == r;
  1475.   int rsgn = x->sgn;
  1476.  
  1477.   long ay = (y < 0)? -y : y;
  1478.   int bw = ay / I_SHIFT;
  1479.   int sw = ay % I_SHIFT;
  1480.  
  1481.   if (y > 0)
  1482.   {
  1483.     int rl = bw + xl + 1;
  1484.     if (xrsame)
  1485.       r = Iresize(r, rl);
  1486.     else
  1487.       r = Icalloc(r, rl);
  1488.  
  1489.     unsigned short* botr = r->s;
  1490.     unsigned short* rs = &(botr[rl - 1]);
  1491.     const unsigned short* botx = (xrsame)? botr : x->s;
  1492.     const unsigned short* xs = &(botx[xl - 1]);
  1493.     unsigned long a = 0;
  1494.     while (xs >= botx)
  1495.     {
  1496.       a = up(a) | ((unsigned long)(*xs--) << sw);
  1497.       *rs-- = extract(down(a));
  1498.     }
  1499.     *rs-- = extract(a);
  1500.     while (rs >= botr)
  1501.       *rs-- = 0;
  1502.   }
  1503.   else
  1504.   {
  1505.     int rl = xl - bw;
  1506.     if (rl < 0)
  1507.       r = Icopy_zero(r);
  1508.     else
  1509.     {
  1510.       if (xrsame)
  1511.         r = Iresize(r, rl);
  1512.       else
  1513.         r = Icalloc(r, rl);
  1514.       int rw = I_SHIFT - sw;
  1515.       unsigned short* rs = r->s;
  1516.       unsigned short* topr = &(rs[rl]);
  1517.       const unsigned short* botx = (xrsame)? rs : x->s;
  1518.       const unsigned short* xs =  &(botx[bw]);
  1519.       const unsigned short* topx = &(botx[xl]);
  1520.       unsigned long a = (unsigned long)(*xs++) >> sw;
  1521.       while (xs < topx)
  1522.       {
  1523.         a |= (unsigned long)(*xs++) << rw;
  1524.         *rs++ = extract(a);
  1525.         a = down(a);
  1526.       }
  1527.       *rs++ = extract(a);
  1528.       if (xrsame) topr = (unsigned short*)topx;
  1529.       while (rs < topr)
  1530.         *rs++ = 0;
  1531.     }
  1532.   }
  1533.   r->sgn = rsgn;
  1534.   Icheck(r);
  1535.   return r;
  1536. }
  1537.  
  1538. IntRep* lshift(const IntRep* x, const IntRep* yy, int negatey, IntRep* r)
  1539. {
  1540.   long y = Itolong(yy);
  1541.   if (negatey)
  1542.     y = -y;
  1543.  
  1544.   return lshift(x, y, r);
  1545. }
  1546.  
  1547. IntRep* bitop(const IntRep* x, const IntRep* y, IntRep* r, char op)
  1548. {
  1549.   nonnil(x);
  1550.   nonnil(y);
  1551.   int xl = x->len;
  1552.   int yl = y->len;
  1553.   int xsgn = x->sgn;
  1554.   int xrsame = x == r;
  1555.   int yrsame = y == r;
  1556.   if (xrsame || yrsame)
  1557.     r = Iresize(r, calc_len(xl, yl, 0));
  1558.   else
  1559.     r = Icalloc(r, calc_len(xl, yl, 0));
  1560.   r->sgn = xsgn;
  1561.   unsigned short* rs = r->s;
  1562.   unsigned short* topr = &(rs[r->len]);
  1563.   const unsigned short* as;
  1564.   const unsigned short* bs;
  1565.   const unsigned short* topb;
  1566.   if (xl >= yl)
  1567.   {
  1568.     as = (xrsame)? rs : x->s;
  1569.     bs = (yrsame)? rs : y->s;
  1570.     topb = &(bs[yl]);
  1571.   }
  1572.   else
  1573.   {
  1574.     bs = (xrsame)? rs : x->s;
  1575.     topb = &(bs[xl]);
  1576.     as = (yrsame)? rs : y->s;
  1577.   }
  1578.  
  1579.   switch (op)
  1580.   {
  1581.   case '&':
  1582.     while (bs < topb) *rs++ = *as++ & *bs++;
  1583.     while (rs < topr) *rs++ = 0;
  1584.     break;
  1585.   case '|':
  1586.     while (bs < topb) *rs++ = *as++ | *bs++;
  1587.     while (rs < topr) *rs++ = *as++;
  1588.     break;
  1589.   case '^':
  1590.     while (bs < topb) *rs++ = *as++ ^ *bs++;
  1591.     while (rs < topr) *rs++ = *as++;
  1592.     break;
  1593.   }
  1594.   Icheck(r);
  1595.   return r;
  1596. }
  1597.  
  1598. IntRep* bitop(const IntRep* x, long y, IntRep* r, char op)
  1599. {
  1600.   nonnil(x);
  1601.   unsigned short tmp[SHORT_PER_LONG];
  1602.   unsigned long u;
  1603.   int newsgn;
  1604.   if (newsgn = (y >= 0))
  1605.     u = y;
  1606.   else
  1607.     u = -y;
  1608.   
  1609.   int l = 0;
  1610.   while (u != 0)
  1611.   {
  1612.     tmp[l++] = extract(u);
  1613.     u = down(u);
  1614.   }
  1615.  
  1616.   int xl = x->len;
  1617.   int yl = l;
  1618.   int xsgn = x->sgn;
  1619.   int xrsame = x == r;
  1620.   if (xrsame)
  1621.     r = Iresize(r, calc_len(xl, yl, 0));
  1622.   else
  1623.     r = Icalloc(r, calc_len(xl, yl, 0));
  1624.   r->sgn = xsgn;
  1625.   unsigned short* rs = r->s;
  1626.   unsigned short* topr = &(rs[r->len]);
  1627.   const unsigned short* as;
  1628.   const unsigned short* bs;
  1629.   const unsigned short* topb;
  1630.   if (xl >= yl)
  1631.   {
  1632.     as = (xrsame)? rs : x->s;
  1633.     bs = tmp;
  1634.     topb = &(bs[yl]);
  1635.   }
  1636.   else
  1637.   {
  1638.     bs = (xrsame)? rs : x->s;
  1639.     topb = &(bs[xl]);
  1640.     as = tmp;
  1641.   }
  1642.  
  1643.   switch (op)
  1644.   {
  1645.   case '&':
  1646.     while (bs < topb) *rs++ = *as++ & *bs++;
  1647.     while (rs < topr) *rs++ = 0;
  1648.     break;
  1649.   case '|':
  1650.     while (bs < topb) *rs++ = *as++ | *bs++;
  1651.     while (rs < topr) *rs++ = *as++;
  1652.     break;
  1653.   case '^':
  1654.     while (bs < topb) *rs++ = *as++ ^ *bs++;
  1655.     while (rs < topr) *rs++ = *as++;
  1656.     break;
  1657.   }
  1658.   Icheck(r);
  1659.   return r;
  1660. }
  1661.  
  1662.  
  1663.  
  1664. IntRep*  compl(const IntRep* src, IntRep* r)
  1665. {
  1666.   nonnil(src);
  1667.   r = Icopy(r, src);
  1668.   unsigned short* s = r->s;
  1669.   unsigned short* top = &(s[r->len - 1]);
  1670.   while (s < top)
  1671.   {
  1672.     unsigned short cmp = ~(*s);
  1673.     *s++ = cmp;
  1674.   }
  1675.   unsigned short a = *s;
  1676.   unsigned short b = 0;
  1677.   while (a != 0)
  1678.   {
  1679.     b <<= 1;
  1680.     if (!(a & 1)) b |= 1;
  1681.     a >>= 1;
  1682.   }
  1683.   *s = b;
  1684.   Icheck(r);
  1685.   return r;
  1686. }
  1687.  
  1688. void (setbit)(Integer& x, long b)
  1689. {
  1690.   if (b >= 0)
  1691.   {
  1692.     int bw = (unsigned long)b / I_SHIFT;
  1693.     int sw = (unsigned long)b % I_SHIFT;
  1694.     int xl = x.rep ? x.rep->len : 0;
  1695.     if (xl <= bw)
  1696.       x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
  1697.     x.rep->s[bw] |= (1 << sw);
  1698.     Icheck(x.rep);
  1699.   }
  1700. }
  1701.  
  1702. void clearbit(Integer& x, long b)
  1703. {
  1704.   if (b >= 0)
  1705.     {
  1706.       if (x.rep == 0)
  1707.     x.rep = &_ZeroRep;
  1708.       else
  1709.     {
  1710.       int bw = (unsigned long)b / I_SHIFT;
  1711.       int sw = (unsigned long)b % I_SHIFT;
  1712.       if (x.rep->len > bw)
  1713.         x.rep->s[bw] &= ~(1 << sw);
  1714.     }
  1715.     Icheck(x.rep);
  1716.   }
  1717. }
  1718.  
  1719. int testbit(const Integer& x, long b)
  1720. {
  1721.   if (x.rep != 0 && b >= 0)
  1722.   {
  1723.     int bw = (unsigned long)b / I_SHIFT;
  1724.     int sw = (unsigned long)b % I_SHIFT;
  1725.     return (bw < x.rep->len && (x.rep->s[bw] & (1 << sw)) != 0);
  1726.   }
  1727.   else
  1728.     return 0;
  1729. }
  1730.  
  1731. // A  version of knuth's algorithm B / ex. 4.5.3.34
  1732. // A better version that doesn't bother shifting all of `t' forthcoming
  1733.  
  1734. IntRep* gcd(const IntRep* x, const IntRep* y)
  1735. {
  1736.   nonnil(x);
  1737.   nonnil(y);
  1738.   int ul = x->len;
  1739.   int vl = y->len;
  1740.   
  1741.   if (vl == 0)
  1742.     return Ialloc(0, x->s, ul, I_POSITIVE, ul);
  1743.   else if (ul == 0)
  1744.     return Ialloc(0, y->s, vl, I_POSITIVE, vl);
  1745.  
  1746.   IntRep* u = Ialloc(0, x->s, ul, I_POSITIVE, ul);
  1747.   IntRep* v = Ialloc(0, y->s, vl, I_POSITIVE, vl);
  1748.  
  1749. // find shift so that both not even
  1750.  
  1751.   long k = 0;
  1752.   int l = (ul <= vl)? ul : vl;
  1753.   int cont = 1;
  1754.   for (int i = 0;  i < l && cont; ++i)
  1755.   {
  1756.     unsigned long a =  (i < ul)? u->s[i] : 0;
  1757.     unsigned long b =  (i < vl)? v->s[i] : 0;
  1758.     for (int j = 0; j < I_SHIFT; ++j)
  1759.     {
  1760.       if ((a | b) & 1)
  1761.       {
  1762.         cont = 0;
  1763.         break;
  1764.       }
  1765.       else
  1766.       {
  1767.         ++k;
  1768.         a >>= 1;
  1769.         b >>= 1;
  1770.       }
  1771.     }
  1772.   }
  1773.  
  1774.   if (k != 0)
  1775.   {
  1776.     u = lshift(u, -k, u);
  1777.     v = lshift(v, -k, v);
  1778.   }
  1779.  
  1780.   IntRep* t;
  1781.   if (u->s[0] & 01)
  1782.     t = Ialloc(0, v->s, v->len, !v->sgn, v->len);
  1783.   else
  1784.     t = Ialloc(0, u->s, u->len, u->sgn, u->len);
  1785.  
  1786.   while (t->len != 0)
  1787.   {
  1788.     long s = 0;                 // shift t until odd
  1789.     cont = 1;
  1790.     int tl = t->len;
  1791.     for (int i = 0; i < tl && cont; ++i)
  1792.     {
  1793.       unsigned long a = t->s[i];
  1794.       for (int j = 0; j < I_SHIFT; ++j)
  1795.       {
  1796.         if (a & 1)
  1797.         {
  1798.           cont = 0;
  1799.           break;
  1800.         }
  1801.         else
  1802.         {
  1803.           ++s;
  1804.           a >>= 1;
  1805.         }
  1806.       }
  1807.     }
  1808.  
  1809.     if (s != 0) t = lshift(t, -s, t);
  1810.  
  1811.     if (t->sgn == I_POSITIVE)
  1812.     {
  1813.       u = Icopy(u, t);
  1814.       t = add(t, 0, v, 1, t);
  1815.     }
  1816.     else
  1817.     {
  1818.       v = Ialloc(v, t->s, t->len, !t->sgn, t->len);
  1819.       t = add(t, 0, u, 0, t);
  1820.     }
  1821.   }
  1822.   if (!STATIC_IntRep(t)) delete t;
  1823.   if (!STATIC_IntRep(v)) delete v;
  1824.   if (k != 0) u = lshift(u, k, u);
  1825.   return u;
  1826. }
  1827.  
  1828.  
  1829.  
  1830. long lg(const IntRep* x)
  1831. {
  1832.   nonnil(x);
  1833.   int xl = x->len;
  1834.   if (xl == 0)
  1835.     return 0;
  1836.  
  1837.   long l = (xl - 1) * I_SHIFT - 1;
  1838.   unsigned short a = x->s[xl-1];
  1839.  
  1840.   while (a != 0)
  1841.   {
  1842.     a = a >> 1;
  1843.     ++l;
  1844.   }
  1845.   return l;
  1846. }
  1847.   
  1848. IntRep* power(const IntRep* x, long y, IntRep* r)
  1849. {
  1850.   nonnil(x);
  1851.   int sgn;
  1852.   if (x->sgn == I_POSITIVE || (!(y & 1)))
  1853.     sgn = I_POSITIVE;
  1854.   else
  1855.     sgn = I_NEGATIVE;
  1856.  
  1857.   int xl = x->len;
  1858.  
  1859.   if (y == 0 || (xl == 1 && x->s[0] == 1))
  1860.     r = Icopy_one(r, sgn);
  1861.   else if (xl == 0 || y < 0)
  1862.     r = Icopy_zero(r);
  1863.   else if (y == 1 || y == -1)
  1864.     r = Icopy(r, x);
  1865.   else
  1866.   {
  1867.     int maxsize = ((lg(x) + 1) * y) / I_SHIFT + 2;     // pre-allocate space
  1868.     IntRep* b = Ialloc(0, x->s, xl, I_POSITIVE, maxsize);
  1869.     b->len = xl;
  1870.     r = Icalloc(r, maxsize);
  1871.     r = Icopy_one(r, I_POSITIVE);
  1872.     for(;;)
  1873.     {
  1874.       if (y & 1)
  1875.         r = multiply(r, b, r);
  1876.       if ((y >>= 1) == 0)
  1877.         break;
  1878.       else
  1879.         b = multiply(b, b, b);
  1880.     }
  1881.     if (!STATIC_IntRep(b)) delete b;
  1882.   }
  1883.   r->sgn = sgn;
  1884.   Icheck(r);
  1885.   return r;
  1886. }
  1887.  
  1888. IntRep* abs(const IntRep* src, IntRep* dest)
  1889. {
  1890.   nonnil(src);
  1891.   if (src != dest)
  1892.     dest = Icopy(dest, src);
  1893.   dest->sgn = I_POSITIVE;
  1894.   return dest;
  1895. }
  1896.  
  1897. IntRep* negate(const IntRep* src, IntRep* dest)
  1898. {
  1899.   nonnil(src);
  1900.   if (src != dest)
  1901.     dest = Icopy(dest, src);
  1902.   if (dest->len != 0) 
  1903.     dest->sgn = !dest->sgn;
  1904.   return dest;
  1905. }
  1906.  
  1907. #if defined(__GNUG__) && !defined(_G_NO_NRV)
  1908.  
  1909. Integer sqrt(const Integer& x) return r(x)
  1910. {
  1911.   int s = sign(x);
  1912.   if (s < 0) x.error("Attempted square root of negative Integer");
  1913.   if (s != 0)
  1914.   {
  1915.     r >>= (lg(x) / 2); // get close
  1916.     Integer q;
  1917.     div(x, r, q);
  1918.     while (q < r)
  1919.     {
  1920.       r += q;
  1921.       r >>= 1;
  1922.       div(x, r, q);
  1923.     }
  1924.   }
  1925.   return;
  1926. }
  1927.  
  1928. Integer lcm(const Integer& x, const Integer& y) return r
  1929. {
  1930.   if (!x.initialized() || !y.initialized())
  1931.     x.error("operation on uninitialized Integer");
  1932.   Integer g;
  1933.   if (sign(x) == 0 || sign(y) == 0)
  1934.     g = 1;
  1935.   else 
  1936.     g = gcd(x, y);
  1937.   div(x, g, r);
  1938.   mul(r, y, r);
  1939. }
  1940.  
  1941. #else 
  1942. Integer sqrt(const Integer& x) 
  1943. {
  1944.   Integer r(x);
  1945.   int s = sign(x);
  1946.   if (s < 0) x.error("Attempted square root of negative Integer");
  1947.   if (s != 0)
  1948.   {
  1949.     r >>= (lg(x) / 2); // get close
  1950.     Integer q;
  1951.     div(x, r, q);
  1952.     while (q < r)
  1953.     {
  1954.       r += q;
  1955.       r >>= 1;
  1956.       div(x, r, q);
  1957.     }
  1958.   }
  1959.   return r;
  1960. }
  1961.  
  1962. Integer lcm(const Integer& x, const Integer& y) 
  1963. {
  1964.   Integer r;
  1965.   if (!x.initialized() || !y.initialized())
  1966.     x.error("operation on uninitialized Integer");
  1967.   Integer g;
  1968.   if (sign(x) == 0 || sign(y) == 0)
  1969.     g = 1;
  1970.   else 
  1971.     g = gcd(x, y);
  1972.   div(x, g, r);
  1973.   mul(r, y, r);
  1974.   return r;
  1975. }
  1976.  
  1977. #endif
  1978.  
  1979.  
  1980.  
  1981. IntRep* atoIntRep(const char* s, int base)
  1982. {
  1983.   int sl = strlen(s);
  1984.   IntRep* r = Icalloc(0, sl * (lg(base) + 1) / I_SHIFT + 1);
  1985.   if (s != 0)
  1986.   {
  1987.     char sgn;
  1988.     while (isspace(*s)) ++s;
  1989.     if (*s == '-')
  1990.     {
  1991.       sgn = I_NEGATIVE;
  1992.       s++;
  1993.     }
  1994.     else if (*s == '+')
  1995.     {
  1996.       sgn = I_POSITIVE;
  1997.       s++;
  1998.     }
  1999.     else
  2000.       sgn = I_POSITIVE;
  2001.     for (;;)
  2002.     {
  2003.       long digit;
  2004.       if (*s >= '0' && *s <= '9') digit = *s - '0';
  2005.       else if (*s >= 'a' && *s <= 'z') digit = *s - 'a' + 10;
  2006.       else if (*s >= 'A' && *s <= 'Z') digit = *s - 'A' + 10;
  2007.       else break;
  2008.       if (digit >= base) break;
  2009.       r = multiply(r, base, r);
  2010.       r = add(r, 0, digit, r);
  2011.       ++s;
  2012.     }
  2013.     r->sgn = sgn;
  2014.   }
  2015.   return r;
  2016. }
  2017.  
  2018.  
  2019.  
  2020. extern AllocRing _libgxx_fmtq;
  2021.  
  2022. char* Itoa(const IntRep* x, int base, int width)
  2023. {
  2024.   int fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
  2025.   char* fmtbase = (char *) _libgxx_fmtq.alloc(fmtlen);
  2026.   char* f = cvtItoa(x, fmtbase, fmtlen, base, 0, width, 0, ' ', 'X', 0);
  2027.   return f;
  2028. }
  2029.  
  2030. ostream& operator << (ostream& s, const Integer& y)
  2031. {
  2032. #ifdef _OLD_STREAMS
  2033.   return s << Itoa(y.rep);
  2034. #else
  2035.   if (s.opfx())
  2036.     {
  2037.       int base = (s.flags() & ios::oct) ? 8 : (s.flags() & ios::hex) ? 16 : 10;
  2038.       int width = s.width();
  2039.       y.printon(s, base, width);
  2040.     }
  2041.   return s;
  2042. #endif
  2043. }
  2044.  
  2045. void Integer::printon(ostream& s, int base /* =10 */, int width /* =0 */) const
  2046. {
  2047.   int align_right = !(s.flags() & ios::left);
  2048.   int showpos = s.flags() & ios::showpos;
  2049.   int showbase = s.flags() & ios::showbase;
  2050.   char fillchar = s.fill();
  2051.   char Xcase = (s.flags() & ios::uppercase)? 'X' : 'x';
  2052.   const IntRep* x = rep;
  2053.   int fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
  2054.   char* fmtbase = new char[fmtlen];
  2055.   char* f = cvtItoa(x, fmtbase, fmtlen, base, showbase, width, align_right,
  2056.             fillchar, Xcase, showpos);
  2057.   s.write(f, fmtlen);
  2058.   delete [] fmtbase;
  2059. }
  2060.  
  2061. char*  cvtItoa(const IntRep* x, char* fmt, int& fmtlen, int base, int showbase,
  2062.               int width, int align_right, char fillchar, char Xcase, 
  2063.               int showpos)
  2064. {
  2065.   char* e = fmt + fmtlen - 1;
  2066.   char* s = e;
  2067.   *--s = 0;
  2068.  
  2069.   if (x->len == 0)
  2070.     *--s = '0';
  2071.   else
  2072.   {
  2073.     IntRep* z = Icopy(0, x);
  2074.  
  2075.     // split division by base into two parts: 
  2076.     // first divide by biggest power of base that fits in an unsigned short,
  2077.     // then use straight signed div/mods from there. 
  2078.  
  2079.     // find power
  2080.     int bpower = 1;
  2081.     unsigned short b = base;
  2082.     unsigned short maxb = I_MAXNUM / base;
  2083.     while (b < maxb)
  2084.     {
  2085.       b *= base;
  2086.       ++bpower;
  2087.     }
  2088.     for(;;)
  2089.     {
  2090.       int rem = unscale(z->s, z->len, b, z->s);
  2091.       Icheck(z);
  2092.       if (z->len == 0)
  2093.       {
  2094.         while (rem != 0)
  2095.         {
  2096.           char ch = rem % base;
  2097.           rem /= base;
  2098.           if (ch >= 10)
  2099.             ch += 'a' - 10;
  2100.           else
  2101.             ch += '0';
  2102.           *--s = ch;
  2103.         }
  2104.     if (!STATIC_IntRep(z)) delete z;
  2105.         break;
  2106.       }
  2107.       else
  2108.       {
  2109.         for (int i = 0; i < bpower; ++i)
  2110.         {
  2111.           char ch = rem % base;
  2112.           rem /= base;
  2113.           if (ch >= 10)
  2114.             ch += 'a' - 10;
  2115.           else
  2116.             ch += '0';
  2117.           *--s = ch;
  2118.         }
  2119.       }
  2120.     }
  2121.   }
  2122.  
  2123.   if (base == 8 && showbase) 
  2124.     *--s = '0';
  2125.   else if (base == 16 && showbase) 
  2126.   { 
  2127.     *--s = Xcase; 
  2128.     *--s = '0'; 
  2129.   }
  2130.   if (x->sgn == I_NEGATIVE) *--s = '-';
  2131.   else if (showpos) *--s = '+';
  2132.   int w = e - s - 1;
  2133.   if (!align_right || w >= width)
  2134.   {
  2135.     while (w++ < width)
  2136.       *--s = fillchar;
  2137.     fmtlen = e - s - 1;
  2138.     return s;
  2139.   }
  2140.   else
  2141.   {
  2142.     char* p = fmt;
  2143.     for (char* t = s; *t != 0; ++t, ++p) *p = *t;
  2144.     while (w++ < width) *p++ = fillchar;
  2145.     *p = 0;
  2146.     fmtlen = p - fmt;
  2147.     return fmt;
  2148.   }
  2149. }
  2150.  
  2151. char* dec(const Integer& x, int width)
  2152. {
  2153.   return Itoa(x, 10, width);
  2154. }
  2155.  
  2156. char* oct(const Integer& x, int width)
  2157. {
  2158.   return Itoa(x, 8, width);
  2159. }
  2160.  
  2161. char* hex(const Integer& x, int width)
  2162. {
  2163.   return Itoa(x, 16, width);
  2164. }
  2165.  
  2166. istream& operator >> (istream& stream, Integer& val)
  2167. {
  2168.   if (!stream.ipfx0())
  2169.     return stream;
  2170.   int sign = ' ';
  2171.   register streambuf* sb = stream.rdbuf();
  2172.   int base = 10;
  2173.   int ndigits = 0;
  2174.   register int ch = sb->sbumpc();
  2175.   while (ch != EOF && isspace(ch))
  2176.     ch = sb->sbumpc();
  2177.   if (ch == '+' || ch == '-')
  2178.     {
  2179.       sign = ch;
  2180.       ch = sb->sbumpc();
  2181.       while (ch != EOF && isspace(ch))
  2182.     ch = sb->sbumpc();
  2183.     }
  2184.   if (ch == EOF) goto eof_fail;
  2185.   if (!(stream.flags() & ios::basefield))
  2186.     {
  2187.       if (ch == '0')
  2188.     {
  2189.       ch = sb->sbumpc();
  2190.       if (ch == EOF) { }
  2191.       else if (ch == 'x' || ch == 'X')
  2192.         {
  2193.           base = 16;
  2194.           ch = sb->sbumpc();
  2195.           if (ch == EOF) goto eof_fail;
  2196.         }
  2197.       else
  2198.         {
  2199.           sb->sputbackc(ch);
  2200.           base = 8;
  2201.           ch = '0';
  2202.         }
  2203.     }
  2204.     }
  2205.   else if ((stream.flags() & ios::basefield) == ios::hex)
  2206.     base = 16;
  2207.   else if ((stream.flags() & ios::basefield) == ios::oct)
  2208.     base = 8;
  2209.  
  2210.   val.rep = Icopy_zero(val.rep);
  2211.  
  2212.   for (;;)
  2213.     {
  2214.       if (ch == EOF)
  2215.     break;
  2216.       int digit;
  2217.       if (ch >= '0' && ch <= '9')
  2218.     digit = ch - '0';
  2219.       else if (ch >= 'A' && ch <= 'F')
  2220.     digit = ch - 'A' + 10;
  2221.       else if (ch >= 'a' && ch <= 'f')
  2222.     digit = ch - 'a' + 10;
  2223.       else
  2224.     digit = 999;
  2225.       if (digit >= base)
  2226.     {
  2227.       sb->sputbackc(ch);
  2228.       if (ndigits == 0)
  2229.         goto fail;
  2230.       else
  2231.         goto done;
  2232.     }
  2233.       ndigits++;
  2234.       switch (base)
  2235.     {
  2236.     case 8:
  2237.       val <<= 3;
  2238.       break;
  2239.     case 16:
  2240.       val <<= 4;
  2241.       break;
  2242.     default:
  2243.       val *= base;
  2244.       break;
  2245.     }
  2246.       val += digit;
  2247.       ch = sb->sbumpc();
  2248.     }
  2249.  fail:
  2250.   stream.set(ios::failbit);
  2251.  done:
  2252.   if (sign == '-')
  2253.     val.negate();
  2254.   return stream;
  2255.  eof_fail:
  2256.   stream.set(ios::failbit|ios::eofbit);
  2257.   return stream;
  2258. }
  2259.  
  2260. int Integer::OK() const
  2261. {
  2262.   if (rep != 0)
  2263.     {
  2264.       int l = rep->len;
  2265.       int s = rep->sgn;
  2266.       int v = l <= rep->sz || STATIC_IntRep(rep);    // length within bounds
  2267.       v &= s == 0 || s == 1;        // legal sign
  2268.       Icheck(rep);                  // and correctly adjusted
  2269.       v &= rep->len == l;
  2270.       v &= rep->sgn == s;
  2271.       if (v)
  2272.       return v;
  2273.     }
  2274.   error("invariant failure");
  2275.   return 0;
  2276. }
  2277.  
  2278. void Integer::error(const char* msg) const
  2279. {
  2280.   (*lib_error_handler)("Integer", msg);
  2281. }
  2282.  
  2283.